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Abstract This paper presents a new method for the reconstruction of weak lensing mass maps. It uses the 
multiscale entropy concept, which is based on wavelets, and the False Discovery Rate (FDR) which allows us 
to derive robust detection levels in wavelet space. We show that this new restoration approach outperforms 
several standard techniques currently used for weak shear mass reconstruction. This method can also be used 
to separate E and B modes in the shear field, and thus test for the presence of residual systematic effects. We 
concentrate on large blind cosmic shear surveys, and illu strate our results using simulated shear maps derived 
from N-Body ACDM simulations llVale a.nd Whit 4 1200,'^! with added noise corresponding to both ground-based 
and space-based observations. 
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1. Introduction 

Weak Gravitational Lensing provides a unique method 
to map directly the distribution of dark matter i n the 
universe j|BartehnamiaiMSclmeideiUl999yMenie^Jl999l 

VanWaerbeke_et_ai [ l200lt iMellieil 120021: iRefregieil 

200,811 . This method is based on the weak distortions that 
lensing induces in the images of background galaxies as 
light travels through intervening structures. This method 
is now widely used to map the mass of clusters and 
superclusters of galaxies and to measure the statistics of 
the cosmic shear field on large scales. 

Ongoing efforts are made to improve the detection 
of cosmic shear on existing telescopes and future instru¬ 
ments dedicated to survey cosmic shear are planned. 
Several methods are used to derive the lensing shear 
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from the shapes of background galaxies. But the shear 
measurements obtained are always noisy, and when it is 
converted into a map of the projected mass k, the result 
is dominated by the noise. 

Several methods have been devised to reconstruct the 
projected mass distribution from the observed shear field. 
The first non-parametric mass re construction was pro- 
Do sed bviKaiser a nd Squill ll 199. i ll a nd further improved 
bviBarte lmannI lll99,^:lKaiseij lll99 ,^ : ISchneider and Seitd 
lll99,’Th : ISoiiires and Kaiseil lll99(1^ . These methods are 
based on linear inversion methods based on smoothing 
with a fixed kernel. Non-linear reconstruction methods 
were proposed usin g a ma xi mum likelihood a p proach 
llSauires and Karsej. 119961 iBartelmann et all Il996t 
|Seitz_et_alJ, [l99^ o r using the max i mum entropy method 
llBridle et al .1 11 9981:iMarshall et allEnfil . 

In this paper, we describe a method for weak lensing 
mass reconstruction based on a wavelet decomposition. 
We use an iterative filtering method with a multiscale 
entropy regularisation to filter the noise. We discuss 
how this decomposition and regularisation functional 
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is particularly well adapted to this problem. In the 
process, we identify significant wavel et coefficient s using 
the False Dis c overy Rate method llMiller et al.L l200lt 
iHookins et al lEH and show how this is superior to 
the standard na thresholding. The FDR method adapts 
its threshold to the features of the data. We concentrate 
on large blind cosmi c shear surveys and u se the ray¬ 
tracing simulations of IVale and Whit3 ll200.lll to test our 
results. We compare the performance of our method to 
Gaussian and Wiener filtering for the reconstruction of 
the mass field in these simulations. We consider conditions 
similar to both ground-based and space-based cosmic 
shear surveys. We also discuss how our method differs 
from other methods based on the maximum entropy prior. 

In section |21 we present the weak shear mass recon¬ 
struction problem. The earlier methods which have been 
proposed to reconstruct the mass map are described in sec¬ 
tion |21 Section 21 presents the Multiscale Entropy method 
and explains why it is a good alternative to standard 
methods. We also propose a modification of the Multiscale 
Entropy, we use the False Discovery Rate (FDR) method 
for detecting the significant wavelet coefficients. A set of 
experiments designed to test our method are described 
section 21 Our conclusions are summarized in section 21 


2. Weak lensing mass reconstruction 

2.1. Weak lensing 


In weak lensing surveys, the shear 7^(0) with i = 
1,2 is derived from the shapes of galaxies at positions 
9 in the image. The shear field 7i(0) can be writ¬ 
ten in terms of the lensing pot ential (see eg. 


Bartelmann and Schneider! lIlQQflTl l 


71 = 

72 = did2'ip, (1) 


where the partial derivatives di are with respect to 9i. 
The convergence k( 0) can also be expressed in terms of 
the lensing potential as 

i^=\ (^1 + dl) ( 2 ) 

and is related to the surface density E(0) projected along 
the line of sight by 


k(0) 


m 


(3) 


where the critical surface density is given by 

^2 D., 


Slcrit — 


47rG DiDis 


(4) 


and G is Newton’s constant, c is the speed of light and Dg, 
Di and Dig are the angular-diameter distances between 
the observer and the galaxies, the observer and the lens, 
and the lens and the galaxies. In practice, the galaxies 
are not at a fixed redshift, and the expressio n for k is an 
averag e of the redshift of the galaxies (see eg. lBartelmannI 
The lensing effect is said to be weak or strong if 
K <C 1 or K > 1, respectively. 

The left panel of Fig.^shows a simulated convergence 
map derived from ray-tracin g through N-body cosm ologi- 
cal simulations performed bv IVale and Whitel 11200,11 . The 
cosmological model is taken to be a concordance ACDM 
model with parameters Dm = 0.3, Da = 0.7, h = 0.7 and 
(7s = 0.8. The simulation contains 512^ particles with a 
box size of 300ft.“^ Mpc. The resulting convergence map 
covers 2x2 degrees with 1024 x 1024 pixels and a assume 
a galaxy redshift of 1. The overdensities correspond to the 
haloes of groups and clusters of galaxies. The rms value of 
K binned in 0.12 arcmin pixels is = 0.023. The typical 
values of k are thus of the order of a few percent, apart 
from the core of massive halos (see figure ^ . The weak 
lensing condition therefore holds in most regions of the 
sky and will be assumed throughout this paper. 


2.2. Mass inversion 


The weak lensing mass inversion problem consists of 
reconstructing the projected (normalized) mass distribu¬ 
tion k{9) from the measured shear field 7i(0) by inverting 
equations o and 121). (Magnification informa tion can also 
be us ed to improve the reconstruction [see iBridle et alJ 
(UmI)], but is typically more noisy than the shear mea¬ 
surements and has not been considered in this paper). For 
this purpose, we take the Fourier transform of these equa¬ 
tions and obtain 


7 i = P^k, 


i = 1,2 


(5) 


where the hat symbol denotes Fourier transforms and we 
have defined = kf + k^ and 


A(k) 

A(k) 


1,2 _ 1.2 
ft/1 h.2 

k^ 

2kik2 


( 6 ) 
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Figure 1. Left: simulated convergence map from (Vale and White, 2003) for a ACDM model. The region shown is 
2x2 square degree. Right: Shear map superimposed on the convergence map , and right shear map. The size and 
direction of each line gives the amplitude and position angle of the shear at this location on the sky. 


with A(^i,fe) = 0 when and A(fci,^ 2 ) = 0 

when ki = 0 or k 2 = 0. 

The shear map 7^ can be calculated from the conver¬ 
gence map K using these expressions. The right panel of 
Fig.Cl shows the shear field associated with the simulated 
convergence field. As is customary, the direction and size 
of the line segment represent the orientation and ampli¬ 
tude of the shear. The rms shear in the 0.12 amin pixels 
of the resulting map is = ( 7 ^ -I- 7|)2 ~ 0.023. 

Note that to recover k from 71 (resp. 72), there is a 
degeneracy when kf = (resp. when fci = 0 or = 0 ). 
To recover k from both 71 and 72, there is a degeneracy 
only when ki = k^ = 0. Therefore, the mean value of k 
cannot be recovered from the shear maps. This is a special 
instance of the well known mass-sheet degeneracy in the 
weak lensing re construction i f only shear information is 
available (see eg lBartelmannl lIlQOfj) for a discussion). 

In practice, the observed shear 7^ is obtained by av¬ 
eraging over a finite number of galaxies and is therefore 
noisy. The relations between the observed data 716,726 
binned in pixels of area A and the true mass map k are 
given by: 

= Pi * K + Ni ( 7 ) 


where Ni and N 2 are noise contributions with zero mean 
and standard deviation (t„ ~ tjj/ ^/Ng, where Ng = UgA is 
the average number of galaxies in a pixel and Ug is the av¬ 
erage number of galaxies per arcmin^. The rms shear dis¬ 
persion per galaxy CTe arises both from measurement errors 
and the intrinsic shape dispersion of galaxies. In this anal¬ 
ysis, we will assume Ue ~ 0.3 as is approximately found 
for ground-based and space-based weak lensing surveys. 
Typical values for the surface density of usable galaxies 
for weak lensing are 

— rig = 20 gal/arcmin^ for ground-based surveys. 

— rig = 100 gal/arcmin^ for space-based surveys. 


From the central limit theorem, this means that for pixels 
with A > 1 amin^, the noise Ni is, to a good approxi- 
mation, Gaussian in both cases and is uncorrelated (see 
iMarshall et a,lJ 1 2002ll for a direct treatment of individual 
galaxy shears using the MEM method). 


2.3. The Inverse Filter: E and B mode 

We can easily derive an estimation of the mass map 
by inverse filtering by noticing that 

A" + A^ = i. (8) 

2(B) 

The least square estimator k; of the convergence k in 
the Fourier domain is: 

= Ablb + P2l2b (9) 

The relation between this estimator and the true mass 

..(E) - ..... 

map is k; = k + N, where N = PiNi + P 2 N 2 . 

Just as any vector field, the shear field 71 (0) can be 
decomposed into a gradient, or electric (E), component, 
and a curl, or magnetic (B), component. Because the weak 
lensing arises from a scalar potential (the Newtonian po¬ 
tential), it can be shown that weak lensing only produces 
B-modes. On the other hand, residual systematics aris¬ 
ing from imperfect correction of the instrumental PSF or 
telescope aberrations, generally generates both E and B 
modes. The presence of B-modes is thus used to test for 
the presence of residual systematic effects in current weak 
lensing surveys. 

The decomposition of the shear field into each of these 
components can be easily performed by noticing that a 
pure B-mode can be transformed into a pure B mode by 
a rotation of the shear by 45°: 71 ^ — 72 , 72 ^ 71 . As a 
result, we can form the following estimator for the B-mode 
“convergence” field 

k[^^ = P27lb -Pi* 72b, (10) 

and check that it is consistent with zero in the absence of 
systematics. 


Figure 2. In the previous region of 2 x 2 square de- 
grees, noisy mass map k) for the same simulation with 
rig = 100 gal/arcmin^, corresponding to space-based ob¬ 
servations. Even in this case, the unfiltered mass map is 
dominated by noise. 
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As follows from equation ||HJ), the noise and 
in K^^^and /tj^^is still Gaussian and uncorrelated. The in- 
verse filtering does not amplify the noise, but kJ mnd 
Itl^^may be dominated by the noise if and 

are large, which is the case in practice. Fig. |21 shows the 
reconstructed mass map using equation when a realistic 
Gaussian noise has been added to the shear maps plotted 
in Fig.QJright. As expected, it is dominated by noise. This 
has motivated the development of different methods in the 
past which we describe below. 


3.2. Maximum Entropy Method 


The Maximum Entropy Method (MEM) is well-known 
and widely used i n image analysis in astronomy 


see 


iBridle et al.1 iIIQq'^: Stardc et al.l ll200lll : iMarshall et alJ 


ll2002ll : IStarck and Murtag-hl ll200^ for a full description). 
It considers both the data and the solution as probability 
density functions and find the solution using a Bayesian 
approach and adding a prior (the entropy) on the solution. 
Several definitions of entr opy exists. The most com mon is 
the definition proposed in lGiiH a.nd Skilling 1 1 99lh : 


3. Earlier Mass Inversion Methods 

3.1. Linear Fiitering 


The standard method I Kaiser and Sauire'3. llQQ.'lll con¬ 
sists in convolving the noisy mass map 
Gaussian window G with standard deviation ac- 


K^^^with 


= G * = G* Pi* 71b + G * P 2 * 726 (11) 


The quality of the resulting estimation depends strongly 

on the value of ao. Fig. |21 shows the variation of the error 

between the original mass map k shown in FiglU and the 
( 

filtered mass map Kq . For this simulation, the optimal 
value of ac lies between 5 and 10 pixels (1 pixel = 0.12 
arcmin) for space observations (i.e. rig = 100 gal/amin^) 
and lies between 20 and 25 pixels for ground observations 
( i.e. rig = 20 gal/amin^). 

An alternative to Gaussian filtering is the Wiener fil¬ 
tering, obtained by assigning the following weight to each 
fc-mode 


w{k) = 


|«(fc)|2 + |iV(fc)P 


( 12 ) 


Where \k{k)\'^ is a model of the true convergence power 
spectrum and is in practice derived from the data. Wiener 
filtering is known to be optimal when both the signal and 
the noise are a realization of a Gaussian Random Field. 
As can be seen from FigHI this assumption is not valid for 
weak lensing mass maps which display non-Gaussian fea¬ 
tures such as galaxy clusters, groups and filaments. Even 
in this case, Wiener Filtering nevertheless leads to rea¬ 
sonable results, generally better than the simple Gaussian 
filtering. 


= X! X! 2/) - m{x, y) - k{x, y) In ( 

m. 11 ' 


Kix,y) 

m{x,y) 


where m is a model, chosen typically to be a sky back¬ 
ground. Hg has a global maximum a,t n = m. MEM does 
not allow negative values in the solution, which is unnatu¬ 
ral for wide field weak shear data or the CMB data, where 
we measure fluctuations around ze ro. To overcome this, it 
has been proposed to replace Hg i Maisinger et 1il l2004ll 
by: 


H+/-{k) = '^^'ip{x,y)-2m- 

'ip{x,y) -f K{x,y) 


X y 

-K{x,y) In 


2 m 


(13) 


where 'ijj{x,y) = \/k^{x, y) + Am?. Here m does not play 
the same role. It is a constant fixed to the expected signal 


rms. 


Mor e generally MEM method presen ts many draw¬ 
backs i Naravan and Nitvananda . llQSfit Starck et all 
l2nnili and various re finements of MEM have been pro¬ 
posed over the yea rs jjWeijJlRQ^ jBontekoeet all Il994l: 
IPantin and Starcil Il99fit IStarck et al.L l200lh . The last 
developm ents have lead to the s o called Mu l tiscale 
Entropy llPantin and StarckL Il996t IStarck et all I 2 OOII: 
IM a,isinger et a,1.L io^J^wEich is based on an undeci- 
mated isotropic wa velet transform (a trous algorithm) 
llStarck et all Il99^ . It has been shown that the main 
MEM drawbacks (model dependent solution, oversmooth¬ 
ing of compact objects, ...) disappear in the wavelet 
framework. A full discussion and comparison between dif¬ 
ferent restoration methods can be found in IStarck et alJ 
( 2 ^. 
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Figures. Reconstruction error as a function of the kernel size ac (in 0.12 amin pixels) for the Gaussian smoothing 
method, with Ug = 20 gal/amin^ (left) and Ug = 100 gal/amin^ (right). 


4. Multiscale Entropy Restoration 

4.1. The Multiscale Entropy 

The Undecimated Isotropic Wavelet Transform 
(UIWT) decomposes an n x n image I as a superposition 
of the form 

,7 

I{k,l) = cjkj + 

1=1 

where cj is a coarse or smooth version of the original 
image I and wj represents the details of I at scale 2~^ 
(see S tarck et al. llStarck et alillQQSHStarck and Miirtaghl 
l2002^ for details). Thus, the algorithm outputs J + 1 
sub-band arrays of size n x n. We will use an indexing 
convention such that j = 1 corresponds to the finest 
scale (high frequencies). T he Multiscale Entropy concept 
jPantin and Star^ llQQdl) consists in replacing the stan¬ 
dard MEM prior (i.e. the Gull and Skilling entropy) by a 
wavelet based prior. The entropy is now defined as 
j-i 

w) = EE h{wj^k.i)- (14) 

1=1 k,i 

In this approach, the information content of an image is 
viewed as sum of information at different scales. The func¬ 
tion h defines the amount of information relative to a given 
wavelet coefficient. Several functions have been proposed 
for h: 

— LOG-MSE: The Multisca le Entropy function used in 
llPantin and StarchIl99(t) (we call it LOG-MSE in the 
following) is defined by: 


T he same multisca l e entr opy function was also derived 
in iMaisinger e~ in il2nn4. 

- NOISE-MSE: In IStarck et all ll200ll) . the entropy is 
derived using a modeling of the noise contained in the 
data: 

r\-wj,k.i\ dh(x) 

= J Pn{\ Wj^k,l I -m)( )x=udu{17) 

where Pn(wj,k,i) is the probability that the coefficient 
Wj^k,i can be due to the noise: Pn{wj,k,i) = Prob(IU >| 
Wj,k,i !)■ For Gaussian noise, we have: 


ey 7 * + 00 

Pn{wj,k,i) = r— / exp(-IU^/2CT|)dIU 
yziraj 

= erfccUaL) 


(18) 


and 


h{wj^k,i) = 




u erfc(- 


U' JO 


V2o 


-)du (19) 


LOG-MSE presents an indeterminacy when the 
wavelet coefficient is equal or close to 0 and the model 
used in equa tion JM. is somewhat ad hoc. This point 
was raised in iMaisinger et all (j20^. A better choice for 
the LOG-MSE wou l d be the Herbert and Leaby function 
I Hebert and Leahvl Il989li (see also the discussion in sec¬ 
tion E|): 


h{wj^k,i) oc log 




( 20 ) 


h{wj^k,i) = - ruj- I Wj^k,i 


'X 


log(^|^)fll5) 


where crx is the total noise standard deviation of the 
data and aj is the noise standard deviation at scale j. 
Km is a user-supplied parameter. 

— ENERGY-MSE: The entropy can be defined as the 
function of the sq uare of the wavelet coefficients 
jStarck et al.Ll200ljl: 


h{'Wj,k,i) 




(16) 


The ENERGY-MSE is quadratic and leads to a strong 
penalization even for wavelet coefficients with high signal- 
to-noise ratio. Such penalization terms are known to over¬ 
smooth the strongest peaks and should not be used for 
the weak leasing mass reconstruction. The NOISE-MSE is 
very close to the li norm (i.e. absolute value of the wavelet 
coefficient) when the coefficient value is large, which is 
known to produce g ood results for the analy sis of piece- 
wise smooth images llPonoho and EladL l200.lll . We there¬ 
fore choose the NOISE-MSE entropy as the most appro¬ 
priate for the weak leasing reconstruction problem. Fig.O 
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shows the li norm penalization function, the ENERGY- 
MSE and NOISE-MSE. The NOISE-MSE penalization 
presents a quadratic behavior for small coefficients and 
a linear one for larger coefficients. More detaifs are given 
in section EH 

4.2. Significant Wavelet Coefficients using the FDR 

In IPantin and Starch dlQQbll , it has been suggested to 
not apply the regularization on wavelet coefficients which 
are clearfy detected (i.e. significant wavefet coefficients). 
The new Multiscale Entropy is: 


hn{wj,k,i) = M{j,k,l)h{wj^k,i) (21) 

where M{j, k, 1) = 1 — M( j, fc, Z) , and M is the multireso¬ 
lution support llMurtag-h et al.L]l995ll : 


M{j,kJ) = 


1 if Wj^k,i is significant 

0 if Wj^k,i is not significant 


( 22 ) 


This describes, in a Booiean way, whether the data con¬ 
tains information at a given scale j and at a given position 
(fc, 1). Commonly, Wj^k,i is said to be significant if the prob¬ 
ability that the wavelet coefficient is due to noise is small, 
i.e. if P{\ W > Wj^k,i I) < e, where P is a given noise distri¬ 
bution function. In the case of Gaussian noise, this amount 
to state that Wj^k,i is significant if | Wj^k,i |> kaj, where aj 
is the noise standard deviation at scale ?, and fc is a con- 
stant, generally taken between 3 and 5 i Murtag-h et al.L 
With this definition, the number of false detections 
depends on both the e value and the image size. 

An alternative approach to this detection strat- 
egy is the False Discovery Rate method (FDR) 
llBeniamini and Hochberel Il99fjl . This technique has re- 
cently be introduced for astronomica l data analysis 
(iMiller et all l200lt iHookins et al It allows us to 

control the ayerage fraction of false detections made oyer 
the total number of detections. It also offers an effective 
way to select an adaptive threshold. The FDR is given by 
the ratio : 


FDR = 


Dn 


(23) 


where Via are the number of pixels truly inactive declared 
active and Da are the number of pixels declared active. 

This procedure controlling the FDR specifies a rate a 
between 0 and 1 and ensures that, on average^ the FDR 


is no bigger than a: 

E{FDR) < y.a < a (24) 

The unknown factor ^ is the proportion of truly inactive 
pixels. A c omplete desc r iption of the FDR method ca n 
be found in iMiller et al.l ll200lll . In lHookins et al.l l|2002ll . 
it has been shown that the FDR outperforms standard 
methods for sources detection. 

Here, we use the FDR method; at each resolution level 
j of the decomposition. We derive a detection thresh¬ 
old Tj (from a value). We have chosen to take a 
different a value per scale. To fix a a value per scale, 
we used The Receiv e r Ope rating Gharacteristic (ROC) 

I Genovese and Fddvl 119911) curves in order to quantify 
the quality of the detection at a given scale for different 
a values. We found that the otj value must increase with 
scale following the relation: Oj = q;o*2 ^ for the spatial ob¬ 
servations and aj = q;o*( 1.7)-^ for the ground observations 
where ao = 0.0125. We then consider a wavelet coefficient. 
Wj^k,i as significant if its absolute value is larger than Tj. 

4.3. Multiscale Entropy Restoration 

Assuming Gaussian noise, the Multiscale Entropy 
restoration method lead to the minimization of the func¬ 
tional. 


J{k) = 


2crS 


( 25 ) 

3 = 1 k,l 


where tT„ the noise standard deviation in k) , J the num¬ 
ber of scales, f3 is the regularization parameter and W is 
the Wavelet Transform operator. The /3 parameter is cal¬ 
culated automatically under the constraint that the resid¬ 
ual should have a standard deviation equal to the noise 
standard deviation. F ull details of the min imization algo¬ 
rithm can be found in IStarck et alJ l|200lll . as well as the 
way to determine automatically the regularization param¬ 
eter (3- 


4.4. Related Work 

4.4.1. The Generalized Wavelet Regularization 

Using a prior suc h that a pixel value is a function of 
its neighborhood Isee lMolina et alJ ll200lh for more details 
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on the Markov Random Field model), the Bayesian solu¬ 
tion consists in adding the following penalization on the 
solution: 


C{k) = - K{x,y + l)f + 

X y 

+ (t>mx,y) - k{x + 1,2/))^)" 


The function cj), called potential function, is an edge pre¬ 
serving function. The term l3J2x^y^(\\ II (^>2/)) 
can also be interpreted as the Gibbs energy of a 
Markov Random Field. Generally, functions (p are cho¬ 
sen with a quadratic pa rt which ensu res a good smooth¬ 
ing of small gradients I Greenl ^9^), and a linear be- 
havior which cancel s the penalization of large gradients 
teouman and Saii^IlQQ,*!!) : 


1. limt^o = Ij smooth faint gradients. 

2. limt^oo = 0, preserve strong gradients. 

3. is strictly decreasing. 

Such functions are often called L 2 -L 1 functions. Examples 
of 4> functions: 


Figure 5. Penalization functions: dashed, li norm (i.e. 
4>{w) =1 w I); dotted I 2 norm (j){w) = ^ (i.e. ENERGY- 
MSE); continuous, multiscale entropy function (NOISE- 
MSE). 


entropy function. We can immediately see that the mul¬ 
tiscale entropy function presents a quadratic behavior for 
small values, and is closer to the h penalization function 
for large values. Penalization function with a I 2 -I 1 behav¬ 
ior are known to be a good choice for image restoration. 

4.4.2. Multiscale MEM and ICF 

_ The multichannel IGF-MEM method llWeirL Il99ll 

Il992h consists in assuming that the visible-space image 
O is formed by a wei^ted sum of the visible-space image 
channels Oj, O = where W is the number of 

channels and Oj is the result of the convolution between 
a hidden image hj with a low-pass filter (IGF) Cj, called 
ICF (Intrinsic Correlation Function) (i.e. Oj = Cj*hj). In 
practice, the ICF is a Gaussian. The MEM-ICF constraint 
is: 


1 . (I>q(x) = x^: quadratic function. 

2. (j)Tv{x) =1 X \: Total Variation. 

3. (p 2 {x) = 2VI + — 2: Hyper-Surface 

(i Ohai ' lnimuer e t al. tr i997 il . 

4. (j> 3 {x) =x‘^/{l + x‘^) (iQeiiiaii and McOluia 

5. ^>4(0:) = 1 — li Peioiia and Malik i l99fl il . 

6 . 4>^{x) = log(l -I- x'^) fli e bw t am - 1 Leah Tii l989 ib 

Figure 0] shows different (p functions. 

It has been shown that this concept can be generalized 
in the wavelet d omain, lead i ng to a multiscale wavelet pe¬ 
nalization term ll.IalobeanuL l200lh : 

Cw{k) = /3 X! *^(ll Up) (26) 

j,k,l 

When (p{x) = x and p = 1, it corresponds to the h norm of 
the wavelet coefficients. In this framework, the multiscale 
entropy deconvolution method is only one special case of 
the wavelet constraint deconvolution method. 

Figure 0 shows the multiscale entropy penalization 
function. The dashed line corresponds to a h penaliza¬ 
tion (i.e. <p(w) =1 w I), the dotted line to a h penaliza¬ 
tion p{w) = and the continuous line to the multiscale 


CicF = '^\hj\ -ruj- I hj I log (27) 

j=i k / 

In iMaisinEfer et al.l ll2004ll . it was argued that the mul¬ 
tiscale entropy is merely a special case of the intrinsic 
correlation function approach, where we replace the ICF 
kernel by a wavelet function. From the strict mathematical 
point of view, this is right, but this vision minimizes com¬ 
pletely the improvement related to the wavelets. All the 
concepts of sparse representation (which is the key of the 
wavelet success in many applications), fast decomposition 
and reconstruction, zero mean coefficients (which allows 
us to get wavelet coefficients which are independent of the 
background and to derive a robust noise modeling) do not 
exist in the ICF-MEM approach. Furthermore, IGF-MEM 
approach requires to estimate accurately the background, 
which may be sometime s a ver y difficult task, and it has be 
shown llBontekoe et al.L Il994ll that the solution depends 
strongly on this estimation. On the contrary. Multiscale 
MEM needs only an estimation of the noise standard de- 
yiation, which is easy to determine. 

For all these reasons, we prefer to keep our vision of 
the multiscale entropy method as a specific case of the 
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Figure 4. Examples of potential function (j). 


generalized wavelet regularization techniques rather than 
as an extension of the ICF approach. 

5. Results 

5.1. Comparison of methods 

We have used a simulated data set obtained using a 
standard A-CDM cosmological model. A part of the k 
mass map and the shear maps is shown in Fig. ^ The 
field size is 2 x 2 square degrees, sampled with 1024 * 1024 
pixels. 

Noisy shear maps, corresponding to both spatial (i.e. 
Ug = 100 gals amin“^) and ground-based observations 
(i.e. Ug = 20 gals amin“^), are created using equation[71 
Then we have reconstructed the two noisy mass maps from 
equation 0 and applied the following methods: 

1. Gaussian filtering with a standard deviation equal to 
ao = 1 amin. 

2. Gaussian filtering with a standard deviation equal to 
(7(5 = 2.5 amin. 

3. Wiener filtering. 

4. Maximum Entropy Method (MEM) using the 
LensEnt2 package. As this code has not been designed 
for manipulating large images, we had to restrict the 
restoration by this method to a field size ofO.5 x 0.5 
square degree, sampled with 256*256 pixels. Since the 
LensEnt2 maps are positivity constrained, as recom¬ 
mended by the author of the LensEnt2 package, we 
have recovered a physical mass by transforming the 
outputs such that the minimum convergence in the 
central quarter of the reconstruction is zero. To op¬ 
timize the IGF, we have maximized the Bayesian evi¬ 
dence value as a function of IGF width, and found that 
maximum evidence is around 210 arcsec for ground 
observations and around 180 arcsec for space observa¬ 
tions. 

5. Multiscale Entropy method. 

The evaluation is done by i) visual inspection of the 
images, ii) calculating the standard deviation between the 
original k mass map and the reconstructed map (i.e. E = 
^"std(k)'^ ), iii) calculating the standard deviation for each 
of their wavelet scales (i.e. Ej = 


iv) calculation the power spectrum of the error E (for 
MEM and multiscale entropy methods). 

The values ac = 1 amin and ac = 2.5 amin have been 
chosen to optimize the Gaussian filtering for rig = 100 gals 
arcmin”^ and rig = 20 gals arcmin“^, respectively. Tabled 
gives the standard deviation of the error for the four recon¬ 
structed mass maps. It shows that i) the Wiener filtering is 
better than the Gaussian filtering and the MEM-LensEnt2 
method and ii) the Multiscale Entropy outperforms the 
three other methods. 

Fig.El shows the error versus the scale (each wavelet 
scale) for both simulations using the Gaussian fil¬ 
tering (continuous line), the Wiener filtering (dotted 
line), the MEM-LensEnt2 filtering (dashed line) and 
the Multiscale Entropy filtering (dotted-dashed line). 
The wavelet scales 1 to 6 correspond to scales of 
0.12,0.23,0.47,0.94,1.87,3.75 amin respectively. We can 
see that the Multiscale Entropy method produces better 
results for all scales. 

Fig.El shows the log power spectrum of the error. It is 
very consistent with the previous one. Indeed, the MEM 
error becomes very important toward the smallest frequen¬ 
cies (largest wavelet scales). The same experiment has 
been done with a smallest IGF (IGF=120 for the spatial 
simulation), but the result is worse, which is not surprising 
since the IGF value was chosen to get the best results. 

Fig.ia shows from top to bottom the reconstructed 
maps for the Gaussian, the Wiener and Multiscale 
Entropy filtering. Fig. (Hlleft corresponds to ground-based 
observations (i.e. rig = 20) and Fig. fright corresponds to 
spatial observations (i.e. rig = 100). 

Fig. El shows the denoising results on a portion of the 
previous image. Fig. El shows the original noise free sim¬ 
ulated image of the 0.5 x 0.5 square degrees field (upper 
left), the Multiscale Entropy Filtering for the spatial sim¬ 
ulated observations {no = 100) (upper right), the MEM- 
LensEnt2 restoration for the ground based observations 
(bottom left) and the spatial observations (bottom right). 

The computation time for the 1024*1024 pixels map is 
4 minutes for the Multiscale Entropy method, 26 seconds 
for the Wiener filtering and 4 seconds for the Gaussian 
smoothing. The computation time for the 256*256 pixels 
map is around 60 minutes (it depends on the convergence 
of the result) using the MEM-LensEnt2 package. 
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Method 

Error {ug = 20 galjamin^) 

Error {rig = 100 galjamin^) 

Gaussian Filtering {aa = 1 amin) 

1.108 

0.775 

Gaussian Filtering [ac = 2.5 amin) 

0.9138 

0.868 

Wiener Filtering 

0.888 

0.770 

MEM-LensEnt2 

1.091 

0.821 

Multiscale Entropy Filtering 

0.888 

0.746 


Table 1. Standard deviation of the reconstruction error with five different methods. 


Figure 6. Standard deviation versus scale for the ground-based simulation (left) and the space-based simulation 
(right). 


Figure 7. Log Power Spectrum of the Error by Multiscale Entropy Filter and MEM for the ground-based simulation 
(left) and the space-based simulation (right). 


5.2. Robustness to missing data 

During the observations, various problem can cause a 
loss of data in the image. For example, it can be due to 
a defect of the camera CCD, generating a dark line or 
a dark row in the image, or to the presence of a bright 
star in the field of view which forces us to remove part 
of the image. In order to study this problem, we mask 
two rectangular areas, setting all pixel values to 0, in the 
shear maps 71 and 72. By inverse filtering, we have derived 
the noisy mass map ki in which we can also visualize the 
lack of data (Fig. EH upper left). Then we have applied 
the three methods, Gaussian filtering, Wiener filtering and 
Multiscale Entropy, to the noisy mass map and the results 
can be seen respectively in Fig. cni upper right, Fig. cni 
bottom left and bottom right. We can see that all three 
methods are robust to the missing data. Note however 
that, for the Wiener filtering, we have assumed perfect 
knowledge of the power spectrum of k, while, in practice, 
its estimation is made more complicated by the complex 
field geometry. 

Fig.ini shows the error versus the scale for both simu¬ 
lations using the Gaussian filtering (continuous line), the 
Wiener filtering (dotted line) and Multiscale Entropy. We 
can see that the Multiscale Entropy still produces better 
results at all scales. Bayesian methods such MEM could 
also take into account properly missing data, however not 
in a straightforward way as when using wavelets. 


5.3. Cluster detection 

Another important aspect of the weak shear mass re¬ 
construction is the possibility to detect clusters and to 
build a catalog. Here, using the FDR in the wavelet space, 
we detect as significant a set of wavelet coefficients. We 
built an isophote map, where each isophote level corre¬ 
sponds to the detection level in a given scale. This isophote 
is overplotted on the true mass map, which allows us to 
visually check the false detections and the missed detec¬ 
tions. A cluster surrounded by two isophotes means that 
it has been detected at two scales. Fig. ^1 left shows the 
isophote map when we use the regular ka thresholding and 
Fig.inileft right shows the isophote map when we use the 
FDR method. We see that the FDR is more sensitive than 
the ka method for the detection, without being contami¬ 
nated by a large number of false detections. Fig.lldlshows 
a zoom of these two maps. Figure fm shows a comparison 
between the Gaussian filtering, the Wiener filtering and 
FDR-Wavelet method for the detection of clusters. In the 
Gaussian and Wiener maps, the isophotes corresponds to 
a ka detection level where fc = 3,4, 5. It shows clearly how 
the FDR-Wavelet method outperforms the other methods. 

5.4. E/B Decomposition 

As explained in E21 a simple diagnostic test for a wide 
range of systematic effects is to search for the presence of 
B-mode in the lensing maps. In order to test it, we have 
simulated mass maps with a B-mode. 
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Figures. Restoration of the 2x2 square degrees ground-based observation (left) and spatial observation (right). 
From top to bottom, Gaussian filtering, Wiener filtering and Multiscale Entropy filtering. 


Figure 9. In a region of 0.5 x 0.5 square degrees, a sixteenth of the original field : Upper left, simulated mass map, 
upper right, Multiscale entropy filtering for Ug = 100 gal/arcmin^. Bottom left, MEM filtering for Ug = 20 gals/amin“^ 
{ICFwidth = 210) and bottom right for Ug = 100 gals/amin“^ {ICFwidth = 180). 


Figure 10. Upper left, noisy shear map (n^ = lOOgal/arcmin^). Upper right, Gaussian filtering. Bottom left, Wiener 
filtering, and bottom right, Multiscale Entropy filtering. 


Figure 16. Noisy simulated mass map 


Fig-Elleft shows a simulated mass map with a lensing 
E-mode signal (left) and an arbitrary B-mode signal (left). 
As usual, we have added a realistic space-based Gaussian 
noise to the shear of this simulation. Fig. M shows the 
noisy mass map resulting. Using the Multiscale Entropy 
filtering, we have then reconstructed the two components 
of the mass map (see ii2.3ll : E-mode in Fig. El left and 
B-mode in Fig. 1171 right. We see clearly that the wavelet 
separation of the E and B modes is very good. Indeed, the 
two main features in the B-mode have well been recovered, 
without interfering with the reconstruction of the E-mode. 

6. Conclusion 

We have presented in this paper a new way to re¬ 
construct weak lensing mass maps. We have modified the 
Multiscale Entropy method in order to take into account 
the FDR. We have shown that this new method outper¬ 
forms several standard techniques currently used for the 
weak shear mass reconstruction. The visual aspect as well 
as objective criteria, such the rms of the error or the rms 
per scale of the error, clearly show the advantages of the 
proposed approach. Experiments have demonstrated that 
it is also robust to missing data. We have also shown that 
a E/B mode separation can also be performed using this 
method, thus providing a useful test for the spatial distri¬ 
bution of residual systematics. Our method allows us also 


to build a catalog of clusters and the use of FDR leads to 
a clear improvement in sensitivity, compared to what has 
been done previously with wavelets. 

Software 

The software related to this paper, MR/Lens, and its 
full documentation are available from the following web 
page: 


\protect\vrule widthOpt\protect\href{http://j starck. 
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Figure 14. The isophotes represent the detected clusters using the Gaussian filtering (upper left), the Wiener hltering 
(upper right) and the wavelet-FDR method. 

Figure 15. left mass map (E-mode), right mass map (B-mode) 

Figure 17. left filtered noisy mass map (E-mode), right filtered noisy mass map (B-mode) 



